
/*******************************************************************************
* Analysis
*******************************************************************************/

***** create tex tables ********************************************************

*** create maintable.tex & maintable_edited

use "$edited/final_data", clear
keep if date==22500

global esttab_opt_nr     "booktabs nomtitles nodepvar nonotes label wrap noomitted collabels(, none) b(a2) se(a2) replace style(tex)" 
global esttab_opt_stats "stats(g_mean s_mean r2_a N, fmt(1 1 2 0) labels("Mean dep. var." "Mean dep. var. sample" "Adjusted R2" "Observations"))"

local munchars "mun_pop mun_pop2 sh_pop_female sh_pop_0_14 sh_pop_15_29 sh_pop_30_44 sh_pop_45_59 sh_pop_60_74 sh_pop_pflichtschule sh_pop_mededu sh_pop_uni salaries pensions sh_pop_foreign district_schwaz"
local attitude "near_vaccination_center near_waldorf sh_fpo"
local prior_covid "sh_cases_pre cum_deaths_27dec2020"
local weight "[iweight=population]"
local pop_res "if gps<=10 & population<=10000"

foreach var of varlist sh_dose_1_20211214 sh_dose_2_20211214 sh_dose_3_20211214 sh_valid_cert_20211214 {
	sum `var' [fweight=population] `pop_res', meanonly
	local s_`var'=r(mean)
	sum `var' [fweight=population], meanonly
	local g_`var'=r(mean)
}

reg sh_dose_1_20211214 sh_antivax_gps `munchars' `attitude' `prior_covid' i.state `pop_res', cluster(municipality_id)

local munchars_lasso "mun_pop sh_pop_0_4-sh_pop_95_99 sh_pop_ahs-sh_pop_foreign sh_ovp-sh_grun salaries pensions sh_pop_foreign turnout i.state i.urban_type_detailed district_schwaz near_vaccination_center near_waldorf"

reg sh_dose_1_20211214 sh_antivax_gps `weight' `pop_res', cluster(municipality_id)
eststo unconditional
estadd scalar s_mean = `s_sh_dose_1_20211214'
estadd scalar g_mean = `g_sh_dose_1_20211214'
local ac=_b[sh_antivax_gps]
local rc=e(r2_a) 

reg sh_dose_1_20211214 sh_antivax_gps `munchars' `weight' `pop_res', cluster(municipality_id)
eststo munchars
estadd scalar s_mean = `s_sh_dose_1_20211214'
estadd scalar g_mean = `g_sh_dose_1_20211214'

reg sh_dose_1_20211214 sh_antivax_gps `munchars' `attitude' `weight' `pop_res', cluster(municipality_id)
eststo munchars_att
estadd scalar s_mean = `s_sh_dose_1_20211214'
estadd scalar g_mean = `g_sh_dose_1_20211214'

reg sh_dose_1_20211214 sh_antivax_gps i.state `weight' `pop_res', cluster(municipality_id)
eststo state
estadd scalar s_mean = `s_sh_dose_1_20211214'
estadd scalar g_mean = `g_sh_dose_1_20211214'

reg sh_dose_1_20211214 sh_antivax_gps `munchars' `attitude' i.state `weight' `pop_res', cluster(municipality_id)
eststo munchars_state
estadd scalar s_mean = `s_sh_dose_1_20211214'
estadd scalar g_mean = `g_sh_dose_1_20211214'
local r=e(r2_a) 

reg sh_dose_1_20211214 sh_antivax_gps `munchars' `attitude' `prior_covid' i.state `weight' `pop_res', cluster(municipality_id)
eststo munchars_state_precov
estadd scalar s_mean = `s_sh_dose_1_20211214'
estadd scalar g_mean = `g_sh_dose_1_20211214'
local a=_b[sh_antivax_gps]
local r=e(r2_a) 

foreach outvar of varlist sh_dose_2_20211214 sh_dose_3_20211214 sh_valid_cert_20211214 {
	reg `outvar' sh_antivax_gps `munchars' `attitude' `prior_covid' i.state `weight' `pop_res', cluster(municipality_id)
	eststo sp_`outvar'
	estadd scalar s_mean = `s_`outvar''
	estadd scalar g_mean = `g_`outvar''	
}

esttab unconditional munchars munchars_att munchars_state munchars_state_precov sp_sh_dose_2_20211214  sp_sh_dose_3_20211214, ///
	se label star(* 0.10 ** 0.05 *** 0.01) drop(_cons) ///
	indicate("Vaccination center = near_vaccination_center" "Demographics = mun_pop mun_pop2 sh_pop_female sh_pop_0_14 sh_pop_15_29 sh_pop_30_44 sh_pop_45_59 sh_pop_60_74" "Education level = sh_pop_pflichtschule sh_pop_mededu sh_pop_uni" "Share foreigners = sh_pop_foreign" "Income = salaries pensions" "District Schwaz = district_schwaz" "Attitudes = near_waldorf sh_fpo" "Prior COVID19 = sh_cases_pre cum_deaths_27dec2020" "State FE = *state*") ///
	mgroups("$1^st$ dose" "$2^nd$ dose" "$3^rd$ dose" "Certificate", pattern(1 0 0 0 0 1 1)) 	///
	noomitted nobase ar2 nomtitles

esttab unconditional munchars_att munchars_state munchars_state_precov sp_sh_dose_2_20211214 sp_sh_dose_3_20211214 ///
    using "$output/maintable.tex", ///
	star(* 0.10 ** 0.05 *** 0.01) ///
	indicate("Municipality characteristics = mun_pop mun_pop2 sh_pop_female sh_pop_0_14 sh_pop_15_29 sh_pop_30_44 sh_pop_45_59 sh_pop_60_74 sh_pop_pflichtschule sh_pop_mededu sh_pop_uni sh_pop_foreign salaries pensions district_schwaz" "State FE = *state*", labels("yes" "")) drop(_cons) ///
	mgroups("$1^{st}$ dose" "$2^{nd}$ dose" "$3^{rd}$ dose", pattern(1 0 0 0 1 1) prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span}))	///
	$esttab_opt_nr $esttab_opt_stats

* adjust formatting for maintable
filefilter "$output/maintable.tex" "$output/maintable_edited.tex", ///
		from("\BStoprule") ///
		to("\BStoprule & \BSmulticolumn{6}{c}{Percent of population with (at least) ...} \BS\BS") replace

		
* Bias correction following Oster(2019)	
local e_bc = `a'-(`ac'-`a')*(1-`r')/(`r'-`rc')	

display "e_bc = `a'-(`ac'-`a')*(1-`r')/(`r'-`rc')"	

display "Bias corrected effect: `e_bc'"		
	
*** create Balance test table **************************************************

* Predict share antivax and homeopathic GPs
global esttab_opt_nr_predict     "booktabs nomtitles nonumber nodepvar nonotes label wide noomitted collabels(, none) b(a2) se(a2) replace style(tex)" 
global esttab_opt_stats_predict "stats(F r2_a N, fmt(2 2 0) labels("F-statistic" "Adjusted R2" "Observations"))"

reg sh_antivax_gps  `munchars' `attitude' ib(9).state `weight' `pop_res', cluster(municipality_id)
eststo prediction_antivax

esttab prediction_antivax using "$output/balance_test.tex", ///
	star(* 0.10 ** 0.05 *** 0.01) ///
	$esttab_opt_nr_predict $esttab_opt_stats_predict
		
		
*** create Robustness table ****************************************************

global esttab_opt_tit    "booktabs nodepvar nonotes label wrap noomitted collabels(, none) b(a2) se(a2) replace style(tex)" 
global esttab_opt_stats_bas "stats(r2_a N, fmt(2 0) labels("Adjusted R2" "Observations"))"

local munchars "sh_cases_pre cum_deaths_27dec2020 near_vaccination_center mun_pop mun_pop2 sh_pop_female sh_pop_0_14 sh_pop_15_29 sh_pop_30_44 sh_pop_45_59 sh_pop_60_74 sh_pop_pflichtschule sh_pop_mededu sh_pop_uni salaries pensions sh_pop_foreign near_waldorf sh_fpo district_schwaz"
local munchars_d "sh_cases_pre cum_deaths_27dec2020 near_vaccination_center mun_pop mun_pop2 sh_pop_female sh_pop_5_9-sh_pop_95_99 sh_pop_pflichtschule sh_pop_mededu sh_pop_uni salaries pensions sh_pop_foreign near_waldorf sh_fpo district_schwaz"

local weight "[iweight=population]"
local pop_res "if gps<=10 & population<=10000"

gen lpop = log(population)

reg sh_dose_1_20211214 sh_antivax_gps `munchars' i.state `weight', cluster(municipality_id)
eststo nopopres

reg sh_dose_1_20211214 sh_antivax_gps `munchars' i.state `weight' `pop_res' & district_schwaz==0, cluster(municipality_id)
eststo no_schwaz

reg sh_dose_1_20211214 sh_antivax_gps sh_brunner `munchars' i.state `weight' `pop_res', cluster(municipality_id)
eststo presidential

reg sh_dose_1_20211214 sh_antivax_gps `munchars_d' i.state `weight' `pop_res', cluster(municipality_id)
eststo munchars_d_state
estadd scalar s_mean = `s_sh_dose_1_20211214'
estadd scalar g_mean = `g_sh_dose_1_20211214'


reg sh_dose_1_20211214 sh_antivax_gps `munchars' i.state `pop_res', cluster(municipality_id)
eststo unweighted

reg sh_dose_1_20211214 sh_antivax_gps `munchars' i.state [iweight=lpop] `pop_res', cluster(municipality_id)
eststo log_weight

reg sh_dose_1_20211214 sh_antivax_gps `munchars' i.state `weight' `pop_res' & gps==1, cluster(municipality_id)
eststo munchars_state_ind

reg sh_dose_1_20211214 sh_antivax_gps `munchars' i.state `weight' `pop_res' & gps>1, cluster(municipality_id)
eststo munchars_state_l1

gen mto_gps = gps>1
reg sh_dose_1_20211214 i.mto_gps##c.sh_antivax_gps `munchars' i.state `weight' `pop_res', cluster(municipality_id)

reg sh_dose_1_20211214 sh_antivax_gps `munchars' i.state `weight' `pop_res' & near_vaccination_center==0, cluster(municipality_id)
eststo center_far

reg sh_dose_1_20211214 sh_antivax_gps `munchars' i.state `weight' `pop_res' & near_vaccination_center==1, cluster(municipality_id)
eststo center_near

reg sh_dose_1_20211214 sh_antivax_docs `munchars' i.state `weight' `pop_res', cluster(municipality_id)
eststo sh_antivax_docs

	
esttab munchars_state nopopres unweighted log_weight munchars_state_ind munchars_state_l1 center_far center_near no_schwaz sh_antivax_docs, ///
	se label star(* 0.10 ** 0.05 *** 0.01) ///
	indicate("State FE = *state*", labels("yes" "")) drop(_cons) ///
	mtitles(Baseline "No pop. res." "Unweighted" "One GP" "> 1 GP" "Center far" "Center near" "Excl. Schwaz" "All physicians") ///
	order(sh_antivax_gps sh_antivax_docs) noomitted nobase ar2
	

esttab munchars_state_precov presidential munchars_d_state nopopres unweighted no_schwaz sh_antivax_docs using "$output/robustness.tex", ///
	star(* 0.10 ** 0.05 *** 0.01) order(sh_antivax_gps sh_antivax_docs) ///
	indicate("Age in 5-year categories = sh_pop_5_9 sh_pop_10_14 sh_pop_15_19 sh_pop_20_24 sh_pop_25_29 sh_pop_30_34 sh_pop_35_39 sh_pop_40_44 sh_pop_45_49 sh_pop_50_54 sh_pop_55_59 sh_pop_60_64 sh_pop_65_69 sh_pop_70_74 sh_pop_75_79 sh_pop_80_84 sh_pop_85_89 sh_pop_90_94 sh_pop_95_99" "State FE = *state*", labels("yes" "")) drop(_cons) ///
	mtitles(Baseline "\shortstack{Antivax \\ voting}" "\shortstack{Detailed \\age}" "\shortstack{No population \\ restriction}" "Unweighted" "\shortstack{Excluding \\ Schwaz}" "\shortstack{All \\doctors}") ///
	$esttab_opt_tit $esttab_opt_stats_bas

filefilter "$output/robustness.tex" "$output/robustness_edit.tex", from("\BSaddlinespace") to("") replace
	

*** create heterogeneity table *************************************************

esttab munchars_state_ind munchars_state_l1 center_near center_far using "$output/heterogeneity.tex", ///
	star(* 0.10 ** 0.05 *** 0.01) order(sh_antivax_gps) ///
	indicate("Municipality characteristics = sh_cases_pre cum_deaths_27dec2020 mun_pop  mun_pop2 near_vaccination_center sh_pop_female sh_pop_0_14 sh_pop_15_29 sh_pop_30_44 sh_pop_45_59 sh_pop_60_74 sh_pop_pflichtschule sh_pop_mededu sh_pop_uni sh_pop_foreign salaries pensions near_waldorf sh_fpo district_schwaz" "State FE = *state*", labels("yes" "")) drop(_cons) ///
	mgroups("Number of GPs" "Vaccination center", pattern(1 0 1 0) prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span}))	///
	mtitles("one GP" "more than one" "near" "far") ///
	$esttab_opt_tit $esttab_opt_stats_bas
	
exit
